#!/usr/bin/perl

##
## Copyright (C) 2009 Luciano Xumerle <luciano.xumerle@gmail.com>
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.
##

use strict;
use warnings;

my %data   = ();
my $reads  = 0;
my $rmulti = 0;
my $len    = 0;
open( FILE, $ARGV[0] ) || die "INSERT VMF FILE!!!\n";
while (<FILE>)
{
    next if ( !m/[RF]\d$/ );
    my $multi = 1;
    $multi = 0 if ( m/\t+([01]):([01]):([01])\t+/ && $1 + $2 + $3 == 1 );
    my @a = split /\t+/;
    $len = length( $a[1] );

    @a = split ',', $a[3];
    foreach my $aaa (@a)
    {
        if ( $aaa =~ m/^(chr.+).fa:(\d+)[RF]/i )
        {
            push @{ $data{$1}->{$2} }, $2 + $len;
        }
    }
    $reads++ if ( $multi == 0 );
    $rmulti++;
}
close FILE;

my %dest  = ();
my $somma = 0;

my $st  = -1;
my $end = -1;
my $ch  = '';
foreach my $chr ( keys %data )
{
    if ( $ch ne $chr )
    {
        $st  = -1;
        $end = -1;
        $ch  = $chr;
    }
    foreach my $start ( sort { $a <=> $b } keys %{ $data{$chr} } )
    {
        my @a = sort { $a <=> $b } @{ $data{$chr}->{$start} };
        if ( $st == -1 )
        {
            $st  = $start;
            $end = $a[$#a];
        }
        else
        {
            if ( $start <= $end )
            {
                $end = $a[$#a];
            }
            else
            {
                print join "\t", ( $ch, $st, $end );
                print "\n";
                $somma += ( $end - $st );
                $st  = $start;
                $end = $a[$#a];
            }
        }
    }
}

my $rapporto = ( $reads * $len ) / $somma;
print "REGIONI COPERTE DA READS (Multi + single): $somma | LUNGHEZZA READ: $len\n";
print "SOMMA READS SINGOLE = $reads | (READS x READLENGTH) / LUNGHEZZA = $rapporto\n";
$rapporto = ( $rmulti * $len ) / $somma;
print "SOMMA READS MULTI = $rmulti | (READS x READLENGTH) / LUNGHEZZA = $rapporto\n";

